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ABSTRACT 

An algorithm is proposed for denoising the signal induced by cosmic strings in the cosmic 
microwave background (CMB). A Bayesian approach is taken, based on modeling the string 
signal in the wavelet domain with generalized Gaussian distributions. Good performance of 
the algorithm is demonstrated by simulated experiments at arcminute resolution under noise 
conditions including primary and secondary CMB anisotropics, as well as instrumental noise. 
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1 INTRODUCTION 

Observations of the cosmic microwave background (CMB) and of 
the Large Scale Structure of the Universe (LSS) have led to the 
definition of a concordance cosmological model. Recently, analy- 
sis of the temperature data of the CMB over the whole celestial 
sphere from the Wilkinson Microwave Anisotropy Probe (WMAP) 
satellite experiment has played a dominant ro le in designi ng this 



preci s e picture of the Unive r se (iBennett et al. 2003; Sperg el et al 



2003: Hinshaw et al 



I2OO7I : ISpergel et all .2007; Hinshaw et al 



2009 ; [Xomatsu et all 20091) . Experiments dedicated to the obser- 



vation of small portions of the celestial sphere have also provided 
their contribution, including the Arcminute Cosmology Bolometer 
Array Receiver (ACBA R) experime nt (Reichar dt et al.ll200 9). the 
Boomerang experiment (I Jones et aL 2006.1 and the Cos mic Back- 
ground Imager (CBI) experiment "(iReadhead et"aDl2004h . 

According to the concordance cosmological model, the cos- 
mic structures and the CMB originate from Gaussian adiabatic 
I perturbations seeded in the early phase of inflation of the Uni- 
verse. However, cosmological scenarios motivated by theories for 
the unification of the fundamental interactions predict the exis- 
tence of topological defects resulting from phase transitions at 
the end of inflation ( Vilenkin & Sh ellard 1994; Hindmarsh 1995a; 
iHindmarsh & Kibblel7l99'5bl : Ixurok & Spergdl ll99Qh . These de- 
fects would have participated to the formation of the cosmic struc- 
tures, also imprinting the CMB. In particular, cosmic strings are a 
line-like version of such defects which are also predicted in th e 
framework of fundamental string theory (Davis & Kibble l2005h . 
As a consequence, the issue of their existence is a central question 
in cosmology today. 



* E-mail: david.hammond@epfl.ch 

t E-mail: yves.wiaux@epfl.ch 

J E-mail: pierre.vandergheynst@epfl.ch 



Cosmic strings are parametrized by a string tension /x, i.e. a 
mass per unit length of string, which sets the overall amplitude of 
the contribution of a string network. Their main signature in the 
CMB is characterized by temperature steps along the string po- 
sitions. Thisjocalized effect, known as the Kaiser-Stebbins (KS) 
effect (E aiser & Stebbins 1984; Gott 1985), hence implies a non- 
Gaussian imprint of the string network in the CMB. The most nu- 
merous strings appear at an angular size around 1 degree on the ce- 
lestial sphere. CMB experiments with an angular resolution much 
below 1 degree are thus required in order to resolve the width of 
cosmic strings. 

Experimental constraints have been obtained on a possible 
string contribution in terms of upper bounds o n the string tension 
/i (Perivolaropoulos 1993; Bevis et al. 2004; Wvman et all 120051 . 
120061 : iBevis et all 120071 : IFraisse et al.ll2007 ). In this context, even 
though observations largely fit with an origin of the cosmic struc- 
tures in terms adiabatic perturbations, room is still available for the 
existence of cosmic strings. 

The purpose of the present work is to develop an effective 
method for mapping the string network potentially imprinted at 
high angular resolution in CMB temperature data, in the perspec- 
tive of forthcoming arcminute resolution experiments. The ob- 
served CMB signal can be modeled as a linear superposition of 
a statistically isotropic but non-Gaussian string signal proportional 
to an unknown string tension, with statistically isotropic Gaussian 
noise comprising the standard component of the CMB induced by 
adiabatic perturbations as well as instrumental noise. 

We take a Bayesian approach to this denoising problem, based 
on statistical models for both the string signal and noise. Our de- 
noising is done in the wavelet domain, using a steerable wavelet 
transform well adapted for representing the strongly oriented fea- 
tures present in the string signal. We show that the string signal co- 
efficients are well described by generalized Gaussian distributions 
(GGD's), which are fit at each wavelet scale using a training simu- 
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lation borrow ed from the set of r ealistic string simulations recently 
produced by Fraisse et 

We develop a Bayesian least squares method for the denoising, 
where each coefficient of the wavelet decomposition is estimated as 
the expectation value of its posterior probability distribution given 
the observed value. This estimation is split into two parts. Assum- 
ing the string tension is known, we use the GGD model to com- 
pute an estimate of the string signal. However, the true string ten- 
sion is unknown. We deal with this by using a preliminary power 
spectral model (PSM) to calculate a posterior probability for the 
string tension. We then construct our overall estimator by weight- 
ing the GGD based estimators using this posterior probability dis- 
tribution for the string tension. Finally, the string network itself can 
be map ped by taking the m agnitude of the gradient of the denoised 
signal dFraisse et alj|2008l) . 

The denoising algorithm that we present may be considered 
as a modular component of a larger data analysis. Firstly notice 
that the PSM might be replaced by any other model allowing com- 
putation of the posterior probability distribution of the string ten- 
sion, notably those which rely on the best experimental bounds 
on the string tension. Secondly, as our method produces a tem- 
perature map of the same size as the input, it may also find use 
as a pre-processing step for other methods for cosmic string de- 
tection based on explicit edge detec tion ( Jeong & Smoot 2005 ; 
ILo & Wrightll2005l : lAmsel et al.ll2007h . 

The performance of our denoising algorithm is evaluated un- 
der different conditions, with astrophysical noise components in- 
cluding various contributions to the standard component of the 
CMB (i.e. primary and secondary CMB anisotropics), as well as 
instrumental noise. Three quantitative measures of performance 
are considered, namely the signal-to-noise ratio, correlation coeffi- 
cient, and kurtosis of the magnitude of gradient of the string signal. 
Our analyses rely on 100 test simulations of a string signal buried 
in the noise. For each string tension and noise condition consid- 
ered, the test simulations are produced b y linear superpositio n of a 
unique test string simulation, also from dFraisse et al.(l2008h . with 
independent noise realizations. 

In each noise condition, we find that the lowest values for the 
string tension down to which our quantitative measures show effec- 
tive denoising are very close to the lowest value where strings begin 
to be visible by eye. Moreover, we acknowledge that this value is 
slightly larger than a detectability threshold set on the basis of the 
PSM. 

The remainder of this paper is organized as follows. In Sec- 
tion |2l we discuss the string signal and the noise at arcminute res- 
olution, and we describe our numerical simulations. In Section [3] 
we describe the steerable wavelet formalism on the plane and study 
the sparsity of the string signal in terms of the wavelet decomposi- 
tion. In Section|4] we define in detail our wavelet domain Bayesian 
denoising (WDBD) algorithm. In Section |5] we evaluate the per- 
formance of the algorithm. We finally conclude in Section [6] 



2 STRING SIGNAL AND NOISE 

In this section, we describe the string signal, discuss current and ex- 
pected future experimental constraints, and detail the various noise 
components at arcminute resolution. We also describe the numeri- 
cal simulations used for our developments. 



2.1 String signal 

In an inflationary cosmological model, the phase transitions re- 
sponsible for the formation of a cosmic string network occur af- 
ter the end of inflation, so as to produce observable defects. From 
the epoch of last scattering until today, the cosmic string net- 
work continuously imprints the CMB. The so-called scaling solu- 
tion for the string network implies that the most numerous strings 
are imprinted just after last scattering and have a typical angu- 
lar size around 1 degree, of the order of the horizon si ze at that 
time f Vachaspati & Vilenkin 1984; Kibble 1985; Albrecht & Turok 
[1985; Bennet 1986; Bennet & Bouchet 1989; Albrecht & Turok 
[1989; Bennet & Bouchet 1990; Allen & Shellard 1990). Longer 
strings are also imprinted in the later stages of the Universe evo- 
lution, but in smaller number, according to the number of corre- 
sponding horizon volumes required to fill the sky. 

The main signature of a cosmic string in the CMB is described 
by the KS effect according to which a temperature step is induced 
in the CMB along the string position. The relative amplitude of this 
step reads as 

^ = (87r7/3) p, (1) 

where /3 — v/c,^ — (1 — is the relativistic gamma factor, 

and p is a dimensionless parameter uniquely associated with the 
string tension /x through 



where G stands for the gravitational constant and c for the speed of 
light. In the following we call p the string tension. 

Analytical models relying on the KS effect and the scaling 
property were defined to simulate the string signal imprinted in the 
CMB. However, in order to produce precise CMB maps account- 
ing for the full non-linear evolution of the string network, one needs 
to resort to numerical simulations. On small angular scales, realis- 
tic simulations can be produced by stacking CMB maps induced in 
different redshift ranges between last scattering and today. The sim- 
ulations we use in this work have been p roduced by this technique 
teouchet etal.lll988l : iFraisse et al.ll2008l) . 

The string signal is understood as a realization of a statisti- 
cally isotropic but non-Gaussian process on the celestial sphere 
with an overall amplitude rescaled by the string tension p, and char- 
acterized by a nearly scale-free angular power spectrum: Cf (/?) = 
p^Ci, where the positive integer index / stands for the angular fre- 
quency index on the sphere. An analytical expression of this spec- 
trum was provided for / larger than a few hundreds by Fraisse et all 
/2OO8), on the basis of their simulations. We consider CMB ex- 
periments with a small field of view corresponding to an angular 
opening r G [0, 27r) on the celestial sphere. In this context, the 
small portion of the celestial sphere accessible is identified to a 
planar patch of size r x r, and we may consider planar signals in 
terms of Cartesian coordinates x = {x,y). The spatial frequen- 
cies may be denoted as A: = (kx^ky) with a radial component 
k = {kl + k'^y^'^. In this Euclidean limit, the radial component 
identifies with the angular frequency on the celestial sphere, be- 
low some band limit set by the resolution of the experiment under 
consideration: I — k < B. Analogously, the planar power spec- 
trum, depending only on the radial component k for a statistically 
isotropic signal, identifies with the angular power spectrum of the 
original signal on the sphere. In particular, for k larger than a few 
hundreds, the nearly scale-free planar power spectrum of the string 
signal s{x) reads as 
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P^k,p) 



2 o-s 



(3) 



withP"(A:) =Cf for/ = /c. 

In this context, the observed CMB signal can be understood as 
a Unear superposition of the string signal and statistically isotropic 
noise of astrophysical and instrumental origin. As will be discussed 
in detail below, this noise is modeled as Gaussian with some an- 
gular power spectrum CJ" . In the Euclidean limit, the correspond- 
ing planar power spectrum for the noise n(x) may be written as 
P^{k) = Cl^ for / = k. The observed noisy signal f{x) is given 
by: 



/ (x) = s (x) + n (x) . 



(4) 



Let us notice that we consider zero mean signals, identifying per- 
turbations around statistical means. 



2.2 Experimental constraints 

Current CMB experiments achieve an angular resolution on the ce- 
lestial sphere of the order of 10 arcminutes, corresponding to a limit 
angular frequency not far above B ~ 10^. At such resolutions, the 
standard component of the CMB primarily contains the Gaussian 
perturbations induced by adiabatic perturbations at last scattering, 
i.e. when the Universe became essentially transparent to radiation. 
These Gaussian anisotropics are referred to as the primary CMB 
anisotropics. In this context, any possible string signal is confined 
to amplitudes largely dominated by these primary anisotropics. The 
constraints mainly come from a best fit analysis of the angular 
power spe ctrum of the overall CM B signal in the WM AP temper- 
ature data JPerivolaroDo ulos 1993 :lBevis et al. 2004; Wvman et al.l 
I2OOI I2OO6I : iBevis et al. .2007.: .Fra isse et al. 2007). The tightest of 



these constraints (Frais se et al.l 12007) gives the following upper 
bound at 68 per cent confidence level: 

P^Pexp = 2.1xlO"^ (5) 

Algorithms have also been designed for the explicit identifica- 
tion of cosmic strings through the observation of the KS effect on 
CMB temperature data. The results of the analysis of the full- sky 
WMAP data typically provide constraints on the string tension two 
orders of magnitude wider than the best fit analy sis of the CMB 
angular power spectrum, i.e. roughly p < 10~^ ('jeong & SmootI 
I2OO5I : |Lo & Wright 2005). The limited angular resolution of the 
WMAP data relative to a typical string width is actually more harm- 
ful for the explicit local detection of cosmic strings than for the es- 
timation of a global parameter such as the string tension through 
the analysis of an angular power spectrum. Corresponding bounds 
have also been studied in the perspective of experiments providing 
high er resolution obser vation of the CMB on small portions of the 
sky (lAmsel et al]l2007h . 



2.3 Noise at arcminute resolution 

Forthcoming experiments will provide access to higher angular 
resolution. The Planck Surveyor satellite experiment will pro- 
vide full-sky C MB data at a r esolution of 5 arcminutes, i.e. with 
B 2x 10 ^ teouched 2004) Fl. The Atacama Cosmology Tele- 
scope (A CT) ( Kosowsky"200^, or the South Pole Telescope (SPT) 
(iRuhl et al. 2004) will map the CMB on small portions of the celes- 
tial sphere at a resolution around 1 arcminute, i.e. with B ~ 10^. 



^ See also Planck Bluebook at h ttp://www.rssd.esa.int/Planckl 
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Figu re 1. Angula r powe r spectra of string signal and noise (borrowed 
from lFraisse et al.l ( l2008h ) as a function of angular frequency in the range 
/ G [10-^,2 X 10"^] in logio ~ ^^g^^Q axes scaling. The spectrum of the 
string signal is represented for a string tension p = 2xl0~^in terms 
of its analytical expression valid at high angular frequencies (red straight 
line). The noise spectra (ordered by decreasing amplitude at low angular 
frequencies) are: the primary CMB anisotropics (black solid line) and the 
gravitational lensing correction (black dashed line), the thermal SZ effect in 
the Ray leigh- Jeans limit (blue solid line), the non-linear kinetic SZ effect 
(green solid line), and the Ostriker-Vishniac effect (magenta solid line). 



A slightly higher resolution might be reached by the radio interfer- 
ometer Arcminute Microkelvin Ima ger (AMI) dJones etal.ll2002l : 
iBarker et al.]l2006l : IZwart et al.ll2008h . We consider the issue of the 
explicit mapping of the string network in the context of such high 
angular resolution experiments. At these resolutions, the so-called 
secondary CMB anisotropics, induced by interaction of CMB pho- 
tons with the evolving universe after last scattering, will dominate 
the primary anisotropics and must be accounted for in the standard 
component of the CMB. 

For the sake of our analyses, we consider that the standard cos- 
mological parameters (i.e. excluding the string tension) are fixed at 
their values in the context of the concordance cosmological model, 
while the string tension remains undetermined. This approxima- 
tion is supported by the already tight experimental bounds ^ on 
the string tension. In other words, we assume that even if the true 
string tension is non-zero it is to be small, and the true values of the 
standard cosmological parameters are close to their present concor- 
dance values. In this context, the angular power spectrum of both 
the primary and secondary anisotropics may be computed on the 
basis of the assumed concordance values for the standard cosmo- 
logical parameters. 

The statistically isotropic Gaussian primary anisotropics ex- 
hibit exponential damping at high angular frequencies. This con- 
trasts with the slow decay of the nearly scale-free angular power 
spectrum of the string signal, which thus dominates over the pri- 
mary anisotropics at high enough angular frequencies. 

The secondary anisotropics include gravitational effects 
such as the Integrated Sachs- Wolfe (ISW) effect, the Rees- 
Sciama (RS) effect, and gravitational lensing, as well as 
re- scattering effects such as the thermal and kinetic Sunyaev- 
Zel'dovich (S Z) effects. The SZ effec ts domi n ate these secondary 



aniso t ropies ( Sunyaev & ZePdovichl I198Q| : iKomatsu & SeliakI 
I2OO2I : iFraisse et al.l l2008h . The ISW and RS effects associated 
with the time evolution of the standard gravitational potentials can 
be neglected at these angular frequencies. One may thus restrict 
the secondary anisotropies considered in the noise to the linear 
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(Ostriker-Vishniac) and non-linear kinetic SZ effects, as well as 
thermal SZ effect. The SZ effects are actually non-Gaussian, spa- 
tially dependent, and the kinetic and thermal effects are correlated. 
As a simplifying assumption, we treat these two effects as two 
independent statistically isotropic Gaussian noise components. 
The effect of gravitational lensing is very small relative to the 
SZ effects, but we still take it into account as a correction to the 
angular power spectrum of the primary anisotropics. 

At arcminute resolution, the thermal and kinetic SZ effects 
have standard deviations around 10 /xK and 5 /xK respectively. They 
also have a slow decay at high angular frequencies and will dom- 
inate the string signal for string tension values below the current 
experimental bound. Arcminute CMB experiments are in fact pri- 
marily dedicated to the detection of these secondary anisotropics. 
Unlike the other effects considered, which have the same black 
body spectrum as the primary anisotropics, the thermal SZ effect 
on the CMB temperature depends on the frequency of observa- 
tion. Its amplitude decreases from the Rayleigh- Jeans limit (null 
frequency) and around 217 GHz where it is expected to vanish, 
before increasing again at higher frequencies. Figure [T] represents 
the angular power spectra as a function of the angular frequency /, 
for a string signal with string tension p = 2xl0~^, the primary 
CMB anisotropics and the correction due to gravitational lensing, 
the Ostriker-Vishniac and non-linear kinetic SZ effect, and the ther- 
mal SZ effect in the Rayleig h- Jeans limi t. These spectra are explic- 
itly borrowed from Fraisse et al.l (l2008 h. again assuming concor- 
dance values for the standard cosmological parameters. 

Instrumental noise also obviously affects signal acquisition. 
Corresponding expected amplitudes for future experiments should 
be lower than the amplitude of secondary anisotropics, b ut still with 
a sta ndard deviation very roughly around 1/iK per pixel (iKosowskvl 
I2OO 6). We will model instrumental noise as Gaussian white noise, 
i.e. with a flat power spectrum. 

In this context, the performance of the denoising algorithm 
to be defined will be studied in the following limits. As a first 
approach, we consider the secondary anisotropics as a statisti- 
cally isotropic Gaussian noise with power spectrum given by the 
Rayleigh- Jeans limit, that is added to the primary anisotropics. One 
can also assume an observation frequency around 217 GHz tak- 
ing advantage of the frequency dependence of the thermal SZ ef- 
fect, and include in the noise secondary anisotropics in absence of 
this effect. This is equivalent to including only the kinetic SZ ef- 
fect and gravitational lensing in the secondary anisotropics. Notice 
that the futur e ACT will have one of its acquisition frequencies 
at 215 GHz jKosowskvl l2006h . In these two cases, instrumental 
noise is considered to be negligible and simply discarded. These 
two different noise conditions are respectively denoted as SA+tSZ 
(secondary anisotropics with thermal SZ effect) and SA— tSZ (sec- 
ondary anisotropics without thermal SZ effect) in the following. 
Analyzing these limits can reveal to what extent the kinetic and 
thermal SZ effect hamper the denoising of the string signal, as a 
function of the string tension. 

Notice that component separation techniques relying on 
the non-Gaussianity of the thermal SZ effect have been de- 
signed for its extraction from the CMB t emperature data, on 



the basis of multi-frequency observations jHobson et al 
Delabroui lle et al.l l2003l : iMaisinger etaDll999l : IPires et al 



1998 



2006 



Bobin et al. 2008). Other component separation techniques relying 



on the non-Gaussianity of the kinetic SZ effect and on its correla- 
tion with the thermal SZ effect have also been proposed for its ex - 
traction from the CMB temperature data (iForni & Aghanimll2004l) . 
In that regard, a global component separation technique might be 



envisaged in order to simultaneously extract all non-Gaussian com- 
ponents of the CMB temperature data, including the string signal. 

In the context of our string signal denoising approach, the per- 
formance of a denoising algorithm can also be examined in the 
limit where the noise only includes primary anisotropics and in- 
strumental noise, assuming secondary anisotropics have been cor- 
rectly separated. The case without instrumental noise is denoted 
as PA— IN (primary anisotropics without instrumental noise) and 
will be studied in order to understand the behaviour of the denois- 
ing algorithm in ideal noise conditions. The case with instrumental 
noise with a standard deviation of IpK, denoted as PA+IN (pri- 
mary anisotropics with instrumental noise), is also considered. 

For the sake of our anal yses, f oreground emissions such as 
Galactic dust or point sources (IKosowskv 2006) are disregarded. 



2.4 Numerical simulations 

Our denoising approach is based on explicitly describing the sta- 
tistical properties of the string signal on small angular scales. We 
need precise simulations of the string signal on planar patches for 
both training and validation of our method. We use 2 simulations 
of the string sig nal borro^yed fr om the full set of 84 simulations 
produced by ( Fr aisse et aDl2008). The first simulation of the string 
signal is used as training data for fitting the prior probability dis- 
tributions for the coefficients of the wavelet decomposition of the 
signal s, while the second is reserved for testing the algorithm. In 
the four noise conditions considered (PA— IN, PA+IN, SA— tSZ, or 
SA+tSZ), this test string signal simulations is combined with 100 
independent realizations of the noise in order to produce multiple 
test simulations. 

The simulations are defined on planar patches of size r x r for 
a field of view defined by an angular opening r = 7.2°. The finite 
size of the patch induces a discretization of the spatial frequencies 
below the band limit B\ k — {2i^/t)p — 50/? with/? = (Px^Py) 
and for integer values px and py with —L^px,Py<L with L = 
rB/27T — B/50. The original maps are sampled on grids with 
2L X 2L uniformly sampled points Xi with 1 ^ i ^ 4L^ for L = 
512. The corresponding pixels thus have an angular size around 
0.42^ The corresponding band limit on kx and ky thus reads B ^ 
2.5x10^. 

The astrophysical and instrumental components of the noise 
are modeled as statistically isotropic Gaussian noise on a planar 
patch with the appropriate power spectra. We consider instrumen- 
tal noise with a standard deviation of l/nK. For each noise com- 
ponent, a simulation may easily be produced by taking the Fourier 
transform of Gaussian white noise, renormalizing each Fourier fre- 
quency value by the square root of the corr esponding power s pec- 
trum, and inverting the Fourier transform (iRocha et al.ll2005h . In 
each noise condition considered, an overall noise simulation is ob- 
tained by simple superposition of the required independent com- 
ponents simulated. The power spectrum of the noise P^{k) is the 
sum of the individual spectra. 

We also include the effect of the experimental beam of a typi- 
cal arcminute experiment in the training string signal simulation, 
as well as in all test simulations for each noise condition con- 
sidered. We simply model this effect by convolution of the string 
signal and astrophysical noise components with a Gaussian kernel 
with a full width at half maximum (FWHM) of 1 arcminute. This 
corresponds to a Gaussian tapering of angular frequencies with a 
FWHM of 2 X 10^, which effectively limits the angular frequencies 
not far above B ~ 10^. Hence the power spectrum P^(k^p) of 
the string signal in relation ^ and the power spectrum of the astro- 
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Figure 2. Simulated maps of the string signal with noise at 1 arcminute resolution, on a field of view of = 1.4°. The top left panel represents the test 
simulation of the string signal. The top middle and top right panels represent the combinations of this string signal simulation for a string tension p = 2x 10~^ 
in the noise conditions PA— IN and SA— tSZ respectively. The bottom panels represent corresponding maps of the magnitude of gradient. 



physical noise components are multiplied by the square modulus of 
the Fourier transform of the experimental beam. The corresponding 
power spectra of the string signal and of the noise in each noise con- 
dition considered are respectively denoted as P*(/c, p) and P^(k). 

For illustration, Figure [3 represents simulated maps of the 
string signal and noise at the resolution considered, as well as cor- 
responding maps of the magnitude of gradient. For visualization 
purposes, we show only one fifth of the total field of view, cor- 
responding to an angular opening r' — 1.4°. At a tension p — 
2x 10~^, with noise only including the primary CMB anisotropics, 
the strings are not visible by eye in the original map itself, while 
part of the network appears in the map of the magnitude of gra- 
dient. This illustrates the natural enhancement of high frequency 
features such as temperature steps by the gradient operator. At the 
same string tension, the presence of the secondary anisotropics 
adds noise at the highest angular frequencies and the strings are 
not visible by eye anymore in either the original map or the map 
of the magnitude of gradient, already when the thermal SZ effect is 
discarded. 
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Figure 3. Steerable Pyramid wavelets with = 6 basis orientations, cen- 
tered at the origin of the plane. From top left to bottom right, the highpass 
filter 7^, the wavelets ^q^j with orientation xn = tt (i.e. q = A^) for the 
spatial scales j = 1 to j = 4, and the lowpass filter ji . 



marginal statistics of the wavelet coefficients of string signal and 
the noise. 



3 WAVELETS AND SIGNAL SPARSITY 

In this section, we firstly describe the steerable wavelet transform 
and reformulate the denoising problem in the wavelet domain. We 
then detail the probability distributions that we use to describe the 



3.1 Steerable wavelets 

Wavelet transforms have become widely used in data analysis and 
image processing in recent years, and have found numerous ap- 
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plications in astrophysics jHobson et al.lll999l : lBarreiro & HobsonI 
l200ll : IStarck et al.l2006ah . In general, a particular transform will be 
useful for signal modeling and denoising if the properties of the sig- 
nal of interest are easier to describe, or more distinct from the noise 
process, in the transform domain than in the original domain. The 
cosmic string signal is characterized by localized, oriented edge- 
like discontinuities. This motivates the use of a transform that is 
well adapted for representing localized oriented features. 

Standard orthogonal wavelets are well localized, but are not 
well suited for arbitrarily oriented features as they have strong bias 
for horizontal and vertical orientations due to their tensor product 
construction. Instead, in this paper we use a steerable wavelet trans- 
form. The transform is parametrized by a number of orientations 
N and spatial scales J. The output of the transform is given by the 
convolution of the original signal / with a set of filters at differ- 
ent scales. These filters are formed by scaling and rotating a single, 
"mother" wavelet 7(x). For the discrete numerical transform, the 
rotation is sampled at N equally spaced angles Xq — Q^/N for 
integer q with 1 ^ ^ ^ A^. The scalings are sampled dyadically, 
i.e. as for integer j with 1 ^ j ^ J. The transform output at 
spatial scale j and orientation q is given by / ★ 7g,j where 7g,j (x) 
is given by rotating 7 by Xg and scaling by 2~^ . 

In order to ensure the invertibility of the transform, it is also 
necessary to include residual highpass and lowpass bands, gener- 
ated by filters 7^ and 7^ respectively. The output of the complete 
steerable wavelet transform then includes a highpass band, J sets 
of N oriented bandpass bands, and the lowpass band. Invertibility 
of the transform is important for our work, as we are interested in 
reconstructing the string map which resides in the image domain. 

We adopt the notation to denote the full vector of wavelet 
coefficients for a given input signal /, where we have implicitly 
vectorized and concatenated the subbands corresponding to differ- 
ent scales and orientations. We will write VK/ to specify individual 
coefficients, where / is a multi-index specifying the scale, orienta- 
tion and spatial location of the coefficient. We will denote by R the 
inverse wavelet transform operator, so that 



(6) 



In this work, we use a particular implementation known as the 
Steerable Pyramicjfl ( Simoncelli et al. 1992). We use the transform 
with = 6 orientations and J = 4 spatial scales. The correspond- 
ing wavelet filters are shown in Figure [3]f or orientation x = tt. In 
particular, the filters that we employ have odd symmetry, which is 
especially appropriate for representing the edge-like discontinuities 
present in the string signal. 

3.2 Problem reformulation in wavelet domain 

By linearity of the wavelet transform, the coefficients of the ob- 
served signal in relation ^ are a sum of the wavelet coefficients 
for the string signal and the Gaussian CMB noise, i.e. 

Wl ^W! + Wi. (7) 

The overall denoising algorithm will proceed by computing the 
wavelet decomposition of the observed signal, estimating the coef- 
ficients corresponding to the string signal, and finally inverting the 
wavelet transform. Our Bayesian estimator requires knowledge of 
probability distributions describing the behaviour of both the string 

See also steerable pyramid implementation available for download at 
|http: //www. ens . nyu . edu/^eero/STEERPYR/, 



signal and the noise. As we shall see later, part of our denoising pro- 
cedure will assume independence of the coefficients for different /, 
allowing it to use a model for the marginal probability distribution 
of the coefficients. 

Notice that by statistical isotropy of both the signal and noise 
processes, the probability distributions of the wavelet coefficients 
for different spatial scales j do not depend on position or orienta- 
tion. For notational convenience we introduce a generalized spatial 
scale b G j, h} for 1 ^ j ^ J. The above comment implies that 
the signal and noise distributions depend only on b. The distribu- 
tion for the string coefficients will depend on the string tension p. 
We write this as the conditional probability Tih{Wf\p). The noise 
coefficient distribution will be denoted gh{WY). 



3.3 String signal distribution 

The morphology of the string signal should give rise to sparse dis- 
tribution of its wavelet coefficients, i.e. many coefficients are close 
to zero with a small number of large magnitude coefficients near the 
temperature steps. We observe this behaviour in our training simu- 
lation. The sparse wavelet coefficients can be successfully modeled 
by a class of probability distributions known as the Generalized 
Gaussian Distributions (GGD). 

We thus use a GGD to model the prior distributions vTb: 



TTb {x\p) 



Vh 



exp 





X 






pub 





(8) 



where F is the Gamma function, and where Vb and pUb are respec- 
tively called shape parameters and scale parameters. These distri- 
butions all have zero statistical means as the signal itself is defined 
in relation (|4l) as a zero mean perturbation. 

Let us acknowledge that GGD's have been used previ- 
ously to model wavelet coefficients for various image process- 
ing app lications inc luding denoising ( Simoncelli & Adel soniri996l : 
Moulin & Liu 1999), deconvolution ( Beig e et al.l200 Q), and coding 
CAntonini et al. 1992; Mallat 1998). 

The shape parameters Vb can be considered as a continu- 
ous measure of the sparsity of the underlying distribution. Setting 
Vb = 2 recovers the Gaussian distribution, which is non-sparse. 
Letting Vb approach yields very peaked probability distributions 
with heavy tails relative to Gaussian distributions, i.e. very sparse 
distributions. These parameters determine the kurtoses nl, i.e. the 
ratio of the fourth central moment to the square of the variance 
(second central moment), by 



F 5^;: 



Kb 



(9) 



The scale parameters pUb are linearly proportional to the stan- 
dard deviations cr^ of the distributions. The corresponding vari- 
ances reflect the power spectrum ^ of the string signal in the range 
of spatial frequencies k probed by the filter at scale 6, and thus also 
scale as p^ : 



{cTb) 



2 
Ub. 



(10) 



The parameters Vb and Ub are estimated by a moment method 
from the wavelet decomposition of the training simulation of the 
string signal for a given string tension p. For each spatial scale, the 
sample variance and kurtosis are calculated, then equations ^ and 
(fTOl ) are solved numerically to obtain Vb and Ub. 
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Figure 4. Logarithm of the modeled prior GGD's tt^ for the wavelet coefficients Wj of a string signal (red solid curves), computed using = 6 orientations 
and J = A spatial scales, superimposed on the histograms of the corresponding coefficients from the training simulation (black dashed curves) for h = {/i, j, 1} 
and 1 ^ j ^ J = 4. From top left to bottom right, the coefficients associated with the highpass filter 7^, with the wavelets 7g,j for the spatial scales j = 1 
to j = 4, and with the lowpass filter ji . 
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cr^ (PA-IN) 


cr^(SA-tSZ) 


cr^(SA+tSZ) 


h 


3.8 xlO^ 


0.41 


1.2 xlO^ 


47 


4.7x10-4 


8.6x10-2 


0.16 


J = 1 


1.6 xlO^ 


0.42 


4.1 xlO^ 


42 


4.5x10-3 


0.26 


0.56 


J = 2 


3.8x10^ 


0.49 


4.5 xlO^ 


26 


0.55 


2.3 


5.4 


J = 3 


7.2 xlO^ 


0.63 


3.1x10'^ 


14 


32 


35 


44 


J = 4 


9.9x10'^ 


0.88 


1.8 xlO^ 


7.3 


4.9x10^ 


4.9 xlO^ 


5.1 xlO^ 




9.8 xlO^ 


1.5 


8.3 xlO^ 


3.7 


2.5 xlO^ 


2.5 xlO^ 


2.5 xlO^ 



Table 1. Parameters for the modeled prior GGD's tt^ and noise Gaussian distributions gij for the spatial scales b = {/i, j, /} and 1 ^ j ^ J = 4. The 
first column identifies the spatial scale b. The next four columns identify the parameters and , and the corresponding standard deviations at p = 1 
and kurtoses nf^ for the prior GGD's. The last three columns identify the standard deviations cr^ for the noise distributions in the noise conditions PA— IN, 
SA— tSZ, and SA+tSZ. All values are given with two significant figures. 



Figure HI shows the modeled prior GGD's 7Vb for the coeffi- 
cients of the string signal with the steerable wavelet 7 with = 6 
orientations and J = 4 spatial scales (see Figure[3]). The GGD's are 
superimposed on the histograms of the corresponding coefficients 
from the training simulation. As the distributions for coefficients of 
different orientations at the same spatial scale will be identical by 
statistical isotropy, the corresponding histograms are produced by 
aggregating the coefficients over all 6 orientations. Qualitatively, 
we see that the prior distributions TTb are well modeled by GGD's, 
which justifies our choice of parameters N and J. 



The estimated values of the parameters Ub and Vb, and cor- 
responding standard deviations and kurtoses k,1 are reported in 
the columns two to five of Table[T] Notice that the shape parameters 
measured for the highpass band (b = h) and for the four bandpass 
bands (j = 1 to j = 4) are significantly lower than 2, correspond- 
ing to very sparse distributions. The larger value for the shape pa- 
rameter for the lowpass band justifies our choice of J = 4 for the 
maximal spatial scale. At the scales accounted for by the lowpass 
filter, the signal coefficients are not significantly non-Gaussian and 
will not be very sparsely distributed. The reconstruction of temper- 
ature steps therefore does not strongly rely on those scales. 



3.4 Noise distribution 

The Gaussian probability distributions gb for the noise coefficients 
W7 are defined as: 



9b (x) 



exp 



(11) 



These distributions are all zero mean, as the noise itself is de- 
fined in relation (lU as a zero mean perturbation. For each of the 
noise conditions PA-IN, PA+IN, SA-tSZ, or SA+tSZ, the vari- 
ances (cr^)^ can be inferred from the power spectrum (k) of the 
noise at 1 arcminute resolution in the range of spatial frequencies 
k probed by the wavelets at the different spatial scales: 



i^b) 



1 

72 ^ 



|7G«rP"(fc). 



(12) 



In this relation, the multi-index value G reads as G = (q^j) with 
1 ^ q ^ N and 1 ^ j ^ J for the oriented wavelet coefficients, 
C = / for the lowpass coefficients, and G = h for the highpass 
coefficients. Notice that due to the rotational invariance of the noise 
power spectrum, the variances calculated in equation f\2\ do not 
depend on the orientation q. The values of the standard deviations 
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cr^ for the noise conditions PA-IN, SA-tSZ, and SA+tSZ are 
listed in the last three columns of Table [T] 



4 BAYESIAN DENOISING 

In this section, we define in detail our wavelet domain Bayesian 
denoising (WDBD) algorithm. We define an overall Bayesian least 
squares estimator as an average of estimation functions evaluated 
at each value of the unknown string tension, weighted by the poste- 
rior probability distribution for the string tension. We then discuss 
Wiener filtering as a standard alternative to our WDBD algorithm. 

4.1 Bayesian least squares 

In a Bayesian approach the signal coefficients are estimated 
from their posterior probability distribution given the coefficients 
of the observed signal : p{W^\W^). Under our general assump- 
tion that the standard cosmological parameters are fixed at their 
concordance values while the string tension remains undetermined 
this posterior probability distribution reads as 

p{w'\W^^ = I dpp{j)\W^^ p{w'\W^,p^, (13) 

where p{p\W^) is the posterior probability distribution function for 
p given and p ( | , p) is the posterior probability distribu- 
tion function for siven and p. 

Several possible methods for selecting an appropriate estimate 
given the posterior probability distribution are possible. Maximiz- 
ing this distribution leads to the maximum a posteriori (MAP) esti- 
mate. Other approaches may consist of minimizing some expected 
cost function. We employ the well known Bayesian least squares 
estimate which minimizes a quadratic cost function. This estimate 
is given by the expectation value of the posterior probability dis- 
tribution. Using relation ( [T3] ) and the linearity of the expectation 
value, our estimator may be written as 

Tr = = j dpp{p\w^^ (14) 

with 

Tr (p) = ^ = I dWW p{w'\w^ ,pY (15) 

This estimation of the signal coefficients is thus given by the mean 
of the estimations for different string tensions, weighted by poste- 
rior probability distribution for p given the observed signal. 

4.2 Estimation functions 

We concentrate firstly on computation of (p) from equation 
iTSl . Notice that the coefficients of the wavelet decomposition of 
the signal and noise are correlated at different orientations, spatial 
scales and positions. Formally one should construct probability dis- 
tributions accounting for these correlations. However, this would 
require computing expectations in a space with dimension equal 
to the number of coefficients of the wavelet decomposition. In this 
perspective, approaches accounting for correlations of the string 
signal devel oped in the framework of maximum en tropy meth- 
ods (MEM) ( Gull & Skilling 1999; Maisinger et al. 1999) might be 
considered. However such methods assume entropic prior models 
and are therefore not directly compatible with our prior model in 
terms of GGD's for the coefficients of the string signal. 



10r 




\wf\ 



Figure 5. Bayesian least squares estimation functions F(-,p) depending 
on observed coefficient Wj at a non-specified spatial scale and for various 
string tensions p. We consider a shape parameter v = 0.5 and various scale 
parameters pu identifying various standard deviations cr* of the coefficients 
of the string signal, all for a unit standard deviation of the noise = 1. 
The black dashed curve shows the limit p ^ oo, where the estimation 
function is the identity function. The upper solid curve (magenta) relates to 
pu = 0.1, i.e. (7^ ~ 1.1, the middle solid curve (blue) relates to pu = 
1.5 X 10~^, i.e. cr^ = 0.16, while the lower solid curve (red) relates to 
pu = 5x 10-3, i.e. = 5.5 x 10-^. 

Accordingly, we employ the simplifying assumption that the 
wavelet coefficients for both the signal and noise for different val- 
ues of / are independent, after conditioning on the string tension p. 
Under this assumption, the integral in expression JTSl l may be fac- 
torized, and each coefficient Wf of the estimate depends only on 
the corresponding observed value Wj . To simplify the notation in 
the following, we write x = Wf and y = Wj to refer to the indi- 
vidual pure string coefficients and observed signal coefficients. We 
shall see that the resulting estimator will also depend on the spatial 
scale b. 

Each coefficient x may then be estimated as a function of the 
corresponding observed value y. By Bayes' theorem we have that 
p{x\y, p) (X p{y\x, p)p{x\p). The probability p{x\p) is exactly the 
marginal probability for each coefficient, which we have modeled 
as the GGD 7Vb. Conditioned on the signal x, the probability of 
observing y is equal to the probability of the noise coefficient being 
equal to exactly y — x. Thus p{y\x, p) is equal to gb{y — x). We 
thus have the posterior probability distribution 

p{x\y, p) = 9b{y - x)'Kb{x\p) (16) 

and our Bayesian estimator at spatial scale b is 

E[x\y, p] = j xgt{y - x)7it{x\p)dx (17) 

with normalization C = f 9b(y — x)7Tb{x\p)dx. This expression 
depends only on the observed coefficient y, the tension p and the 
scale b. This defines the estimation function Fb{y, p) = E[x\y, p]. 
Returning to our original notation, the estimated string coefficients 
are given by evaluation of the estimation function at each scale, i.e. 

W!{p)^Ft(wl,py (18) 

In practice, these estimation functions are computed by nu- 
merical integration and tabulated for the different spatial scales b 
and the required range of string tensions. Figure [5] shows generic 
shapes of estimation functions F{-,p) at a non- specified spatial 
scale and for various string tensions p. For the sake of illustration, 
we consider a shape parameter ^; = 0.5 and various scale param- 
eters pu identifying various standard deviations a* of the coeffi- 
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cients of the string signal, all for a unit standard deviation of the 
noise cr^ = 1. Notice that the estimation functions are odd, and al- 
ways shrink the magnitude of their input, i.e. p) \ < \Wj\. 
Qualitatively, they behave as a smooth thresholding operation on 
the observed coefficient sending small magnitude coeffi- 
cients closer to zero while preserving large magnitude coefficients. 
For small string tensions, the noise dominates the signal and the ef- 
fective thresholding is more severe, while for large string tensions 
the noise becomes negligible and Fb{-, p) reduces to the identity. 

In the particular case of Gaussian signal coefficients (v = 
2), the Bayesian least squares estimation is equivalent to simple 
Wiener filtering of the coefficients. Also notice that in the case of 
Laplacian signal coefficients (v — 1), the estimation function for 
MAP estima tion would reduce to the well known soft-thresholding 
operation ( Moulin & Liull99 ^. By definition, this specific instance 
of thresholding operation sends to zero coefficients with an abso- 
lute value below some threshold, and reduces the absolute value of 
coefficients above the threshold by the value of the threshold itself. 



4.3 Posterior string tension distribution 

By B ayes' theorem, the posterior probability distribution function 
for p given the observed signal p{p\W^) is simply obtained from 
the likelihood C{W^ \p) and the prior probability distribution func- 
tion p{p) on p. For complete consistency, the likelihood should 
be calculated using the framework of the model established for 
the coefficients, based on relation ([7]) and on the prior GGD's vTb. 
However, while this model by construction accounts for the non- 
Gaussianity, i.e. sparsity, of the string signal, it ignores the correla- 
tion between coefficients. 

We have observed that a likelihood yielding a more precise 
localization of the string tension value can actually be obtained us- 
ing a power spectral model. Such a model assumes both the string 
signal and the noise arise from statistically isotropic Gaussian ran- 
dom processes, such that their Fourier coefficients are indepen- 
dent Gaussian variables. This assumption relies on the idea that the 
characteristic temperature steps of the string signal are smoothed 
by projection on the non-local imaginary exponentials defining the 
Fourier basis. Under this model, as the string signal and noise are 
independent, the observed signal / has a power spectrum 



P{k,p)^p- {k) + P^ {k,p), 



(19) 



In this setting the likelihood can be computed most easily in 
terms of the Fourier transform / of the observed signal. Accounting 
for the complex value of the Fourier coefficients as well as for the 
symmetry f{—k) = f*{k) that holds for real signals /(x), this 
likelihood reads as: 



n 

{Px,Py}^ 



1 



\/^P{Kp) 



exp 



1 l/WI 



2p{k,p) 



(20) 



where | • | stands for the modulus of a complex variable. The poste- 
rior probability distribution function for p given the observed signal 
thus reads as 



^[p\f)=D-'p{p)c[f\p) 



(21) 



with normalization D — J p(p)C(f\p)dp. We take the prior p{p) 
to be flat in an interval p G [0, pmax] , with an upper bound pmax large 
enough relative to the upper bound associated with the best exper- 
imental constraints ([5]): pmax > Pexp- In practice, C{ f\p) decays so 
rapidly for large p that the resulting posterior is not sensitive to the 



value pmax provided that it is greater than the effective support of 

^(/Ip). 

We use this PSM posterior p{p\ f) in the place of p{p\W-^) in 
equation ilAi . Each component of the string coefficient is thus 
estimated as 



wf = Jp{p\m{wf., 



p)dp. 



(22) 



In practice, this integral is computed numerically by sampling 20 
values of p chosen to cover the effective support of p{p\f). 

The estimated string signal in the original image domain is 
then given by inverting the wavelet transform, i.e. 



s(x) - [Ryf] {x) 



4.4 Alternative Wiener filtering 



(23) 



In order to obtain a more precise estimation of the posterior prob- 
ability distribution function for p, we have explicitly set up a PSM 
assuming that the string signal arises from a statistically isotropic 
Gaussian random process, such that its Fourier coefficients are in- 
dependent Gaussian variables, just as for the noise. 

At each string tension allowed by p(p|/), one may now con- 
sider estimating the string signal s from the observed signal / sim- 
ply using this Gaussianity assumption. In this case the Bayesian 
least squares estimate for a string tension p reduces to Wiener fil- 
tering in the Fourier domain, so that: 



s{k,p) 



(24) 



[k) + P' {k.pY 
Analogously to relation ([22]), the estimate of the string signal 



in the Fourier domain is 

s{k) = j p{p\f)s{k,p)dp, 



(25) 



and the estimate in the image domain is recovered by inverting the 
Fourier transform. 

Let us acknowledge the fact that this alternative Wiener fil- 
tering based procedure relies only on the knowledge of the power 
spectra of both the signal and noise, while our WDBD approach 
relies on a training simulation for an explicit modeling of the 
prior GGD's for the coefficients of a wavelet decomposition of 
the string signal. However, from the theoretical point of view it is 
clear that the Wiener filtering approach, which disregards the non- 
Gaussianity of the signal to be recovered, will be less effective at 
identifying this signal than our WDBD procedure, which explicitly 
accounts for the corresponding sparsity. While the Gaussianity as- 
sumption is useful for estimating a single global parameter such as 
the string tension on the basis of a PSM, it is not optimal for the 
explicit reconstruction of the sparse features of the string network 
itself. This fact is illustrated in our analysis the algorithm perfor- 
mance in the next section. 



5 ALGORITHM PERFORMANCE 

In this section, we firstly define the WDBD performance criteria 
to be the signal-to-noise ratio, correlation coefficient, and kurtosis 
of the map of the magnitude of gradient of the string signal. We 
then study the algorithm performance in each noise condition, in 
comparison with Wiener filtering. We also examine a detectability 
threshold on the string tension based on on the PSM, and compare 
it with an eye visibility threshold for the WDBD algorithm. 
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5.1 WDBD performance criteria 

As already emphasized, denoising may be used as a pre-processing 
step for other methods for cosmic string de t ection based on ex- 
pUcit edge detect ion (Jeong & SmooJ l2005l : ILo & Wrighl l2005l : 
lAmsel et al The relative performance of such methods be- 

fore and after denoising might be an effective criterion for evaluat- 
ing the denoising performance itself. Here we evaluate the perfor- 
mance of WDBD independently of any further processing. 

The overall denoising is effective for string mapping if the 
magnitude of gradient of the denoised signal closely resembles the 
magnitude of gradient of the true string signal. A simple qualita- 
tive measure of the denoising performance is given by whether the 
string network is visible in the magnitude of gradient of the de- 
noised signal. We define the eye visibility threshold as the mini- 
mum string tension around which the overall denoising and map- 
ping by the magnitude of the gradient begins to exhibit string 
features visible by eye. We will augment this qualitative assess- 
ment of the denoising performance with three quantitative mea- 
sures, namely the signal-to-noise ratio, the correlation coefficient 
and the kurtosis of the magnitude of the gradient of the denoised 
string signal. The first two of these are computed with respect to 
the original known signal, while the kurtosis is computed only us- 
ing the denoised signal. The kurtosis is known to be a good statistic 
for discriminating between models with and without cosmic strings 
(^oessner et al. 1994). 

The signal-to-noise ratio is defined in terms of the magnitude 
of gradient | Vs| (x) of the original string signal s(x) in relation ©, 
and of the magnitude of gradient |Vs|(x) of the denoised signal 

in relation (l23l l as 



SNR' 



(|Vs|,|Vs|) 



-20 log 1 



^(|Vs|-|Vs|) 



(26) 



where cr*^!^*'"'^*!^ and cr'^*' respectively stand for the standard 
deviations of the discrepancy signal |Vs| — |Vs| and of the orig- 
inal signal |Vs|. The standard deviations are estimated from the 
sample variances on the basis of the signal realizations concerned. 
With this definition, the SNR^'^^I'I^^I^ G M is measured in deci- 
bels (dB). Large negative and positive values are respectively asso- 
ciated with large and small discrepancy signals relative to the orig- 
inal signal. An exact recovery of the string network would provide 
an infinite signal-to-noise ratio. We will consider that the denois- 
ing is effective in terms of signal-to-noise ratio for the values of p 
where this statistic is larger after denoising than before, and posi- 
tive. 

The correlation coefficient is defined in terms of the magni- 
tude of gradient of the original and denoised string signals as 



J|V5|,|V^|) 



:ov(l^-l'l^^l) 



(27) 



where cov^'^*' '^^'^ stands for the covariance between |Vs| and 
I Vs|. This signal covariance is also estimated from the sample co- 
variance on the basis of the signal realizations concerned. An exact 
recovery of the string network would provide a unit correlation co- 
efficient. The null value corresponds to a reconstruction completely 
decorrelated from the original signal. We will consider that the de- 
noising is effective in terms of correlation coefficient for the values 
of p where this statistic is larger after denoising than before, and 
positive. 

Analogously, kurtoses are estimated from the sample kurtoses 
on the basis of the signal realizations concerned. The estimated kur- 
tosis of the magnitude of gradient of pure Gaussian noise is dis- 



tributed around a mean value across all test simulations a^'^^' ~ 3, 
even though the magnitude of gradient itself is not Gaussian. At the 
arcminute resolution considered, the estimated kurtosis of the mag- 
nitude of gradient of a pure string signal is much higher than the 
value associated with pure noise, with a value around k'^^' ~ 32 in 
the training simulation. The estimated kurtosis of the magnitude of 
gradient of a string signal with noise before denoising naturally lies 
in the interval [z^:'^^' , a^'^^'], for any value of the string tension p. 
Its mean value a^'^-^' (p) obviously increases fromz^l^^l forp = 
to At'^*' for p ^ oo. An ideal denoising procedure should recover 
exactly the original string signal. The mean value of the estimated 
kurtosis would then be raised to a^'^*' after denoising. In practice, 
the estimated kurtosis of the magnitude of gradient after denois- 
ing is distributed around some mean value n\^'\p) as a function 
the string tension p. The comparison a^'^*' {p) after denoising with 
(p) before denoising measures the denoising performance as 
a function of the string tension. We will simply consider that the 
denoising is effective in terms of kurtosis for the values of p where 
this statistic is significantly larger after denoising than before. 

Our denoising experiments for each noise condition consid- 
ered are performed for string tensions equi- spaced in logarithmic 
scaling in the range logio P ^ ["lO? —05], corresponding to ratio 
values for p of 1.0, 1.6, 2.5, 4.0, and 6.3 in each order of mag- 
nitude. For each noise condition and string tension considered, we 
perform 100 denoising simulations at 1 arcminute resolution. We 
consider that the quantitative measures described above indicate ef- 
fective denoising performance for a given string tension when they 
show effective performance significant over the entire ensemble of 
denoising simulations. 

5.2 Noise conditions PA-IN and PA+IN 

For the PA— IN condition, the magnitude of gradient of the string 
signal before and after WDBD and Wiener filtering is represented 
in Figure [6] for various string tensions from a single simulation. 
Only one fifth of the total field of view of the simulations is shown, 
corresponding to an angular opening r' — 1.4°. 

The visibility of the individual strings of the network is clearly 
enhanced by the denoising. For a value of the string tension around 
the experimental upper bound y9 = 2.5xl0~^, part of the network 
is visible by eye before denoising. At tension p — 6.3 x 10~^, a 
very reduced number of strings is visible by eye before denoising. 
For both of these string tensions, part of the network is visible by 
eye through WDBD and Wiener filtering, but the resulting map is 
clearly more noisy in the second case. The value y9 = 6.3xl0~^° 
is the lower bound on the string tension where a very reduced num- 
ber of strings is visible by eye through WDBD, while no strings 
are visible by eye before denoising. In this limit, only string loops 
are actually recovered, together with some spurious point sources. 
Wiener filtering only provides noise at that level. 

The posterior probability distributions for the string tension 
are reported in Figure |7] as computed from the signals observed at 
the three string tensions of interest in Figure[6] The graphs highlight 
the high precision of the localization of p by the PSM described in 
Section |4] The slight offset observed is not related to a bias of the 
procedure itself but is simply due to an effective difference between 
the power spectrum of the test string simulation and the analytical 
expression of the power spectrum P^{k) used in relations (|20l ) and 
(l2n ). This difference may be associated with a cosmic variance in- 
cluding the contribution of a string signal. 

The signal-to-noise ratio, correlation coefficient, and kurtosis 
of the magnitude of gradient of the string before and after WDBD 
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Figure 6. Magnitude of the gradient of the string signal before denoising (left panels), after Wiener filtering (middle panels) and WDBD (right panels), in the 
noise conditions PA— IN and for various string tensions at 1 arcminute resolution on a field of view of = 1.4°. From top to bottom, the string tensions 
considered are p = 2.5 x 10"'^, p = 6.3 x IQ-^, and p = 6.3 x 10-1°. 



and Wiener filtering are represented in Figure [8] as functions of 
the string tension. In the range of string tensions where the de- 
noising procedure provides visibility of strings by eye, it appears 
clearly that WDBD and Wiener filtering both significantly increase 
the signal-to-noise ratio and correlation coefficient to strictly posi- 
tive values. At low string tensions the correlation coefficient is also 
significantly higher for WDBD than for Wiener filtering. This rep- 
resents a first quantitative measure of the superiority of our ap- 
proach. The kurtosis of the magnitude of gradient is also signif- 
icantly increased from its value before denoising towards higher 
values through WDBD. The peak obtained at low string tensions, 
with kurtosis values above the expected value around /t:'^*' 32, 



reflects the fact that the denoising recovers a thresholded version 
of the string signal in that range, only keeping localized loops in 
the limit identified by the visibility by eye (see Figure [6]). At low 
string tensions, Wiener filtering essentially fails to increase the kur- 
tosis values towards the expected value. We interpret this failure as 
a quantitative measure of the fact that Wiener filtering fails to re- 
move a substantial part of the noise, in contrast with WDBD. This 
represents a second quantitative measure of the superiority of our 
approach. Let us emphasize that the lowest string tensions where 
each of our quantitative measures begin to show effective denois- 
ing performance for the WDBD algorithm are very close to the eye 
visibility threshold. 
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Figure 7. Posterior probability distributions (red solid curves) for the string tension as computed from simulated signals observed for p = 2.5 x 10~^ (left 
panel), p = 6.3xl0~^ (middle panel), and p = 6.3 x 10~^^ (right panel), in the noise conditions PA— IN at 1 arcminute resolution. The black dashed 
vertical lines represent the exact values of the string tension relative to the test string simulation. 
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Figure 8. Signal-to-noise ratio (left panel) in decibels (dB), correlation coefficient (middle panel), and kurtosis (right panel) of the magnitude of gradient 
as functions of the string tension in logarithmic scaling in the range logio P ^ [~^^^ —05] in the noise conditions PA— IN and at 1 arcminute resolution. 
The black dashed curves represent values before denoising, while the red solid curves and green dot-dashed curves represent values after WDBD and Wiener 
filtering respectively. The vertical lines on the curves represent the variability at one standard deviation of the estimated statistic across the 100 test simulations 
considered (these lines are not visible where smaller than the width of the curves). The blue dotted vertical lines represent the eye visibility threshold 
p = 6.3xl0~^^. The blue dotted horizontal lines identify either the limit of zero signal-to-ratio, zero correlation coefficient, or the kurtosis of the magnitude 
of gradient of a pure string signal: ^cl^^l ~ 32. 



The degradation of the denoising performance due to instru- 
mental noise is probed in the noise condition PA+IN, with an in- 
strumental noise level of l/j^K. For this case we omit a complete 
analysis of all our quantitative measures. We simply notice that 
such a small level of instrumental noise already significantly affects 
the denoising performance by raising the eye visibility threshold 
by more than one order of magnitude. The only reason why an ef- 
fective reconstruction of strings may be achieved down to so small 
string tensions in the noise conditions PA— IN is simply that, at high 
spatial frequencies, the string signal with a nearly scale-free power 
spectrum largely dominates the primary anisotropics with an expo- 
nentially damped power spectrum. This advantage is lost as soon 
as high frequency noise is added, in particular instrumental noise. 



5.3 Noise conditions SA-tSZ and SA+tSZ 

The magnitude of gradient of the string signal before and after 
WDBD and Wiener filtering is represented in Figure [9] for various 
string tensions, from a single simulation. In the noise conditions 
SA— tSZ and for a value of the string tension around the experimen- 
tal upper bound p — 4.0x10"^, a very reduced number of strings is 
visible by eye before denoising. Part of the network is visible by eye 
after WDBD and Wiener filtering, but the resulting map is clearly 
more noisy in the second case. The value p — 1.0 x 10~^ is the 
lower bound on the string tension where a very reduced number of 
strings is visible by eye through WDBD, while no string is visible 



by eye before denoising. Wiener filtering only provides noise at that 
level. In the noise conditions SA+tSZ, the value p — 2.5xl0~^is 
the lower bound on the string tension where a very reduced num- 
ber of strings is visible by eye through WDBD, while no string is 
visible by eye before denoising. Again, Wiener filtering only pro- 
vides noise at that level. At the lower bounds for the string tensions 
in both noise conditions only string loops are recovered, still with 
some spurious point sources. 

The posterior probability distributions for the string tension 
are reported in Figure [TOl as computed from the signals observed in 
the three cases of interest in Figure [9l The graphs still highlight the 
high precision of the localization of p by the PSM. 

The signal-to-noise ratio, correlation coefficient, and kurtosis 
of the magnitude of gradient of the string signal before and after 
WDBD and Wiener filtering are represented in Figure [TT] as func- 
tions of the string tension. As for the PA— IN and PA+IN cases, 
both WDBD and Wiener filtering increase the signal-to-noise ra- 
tio and correlation coefficient to strictly positive values for ten- 
sions above the eye visibility threshold. The correlation coefficient 
is significantly higher for WDBD than for Wiener filtering in the 
whole range of string tensions of interest. As before, the kurtosis 
of the magnitude of gradient is also significantly increased from 
its value before denoising towards higher values through WDBD, 
with a peak at low string tensions due to the fact that the denoising 
recovers a thresholded version of the string signal. Wiener filtering 
essentially fails to increase the kurtosis values towards the expected 
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Figure 9. Magnitude of the gradient of the string signal before denoising (left panels), after Wiener filtering (middle panels) and WDBD (right panels), for 
various string tensions at 1 arcminute resolution on a field of view of = 1.4°. The two top panel rows relate to the noise conditions SA— tSZ for string 
tensions p = 4.0x10"^ and p = 1.0x10"^ respectively. The bottom panel row relates to the noise conditions SA+tSZ for a string tension p = 2.5xl0~^. 



value in the whole range of string tensions of interest, once more 
reflecting its poorer denoising performance. As before, we see that 
for both SA— tSZ and SA+tSZ, the eye visibility thresholds are 
very close to the lowest string tensions where each of our quantita- 
tive measures begin to show effective denoising performance. 

Let us acknowledge the fact that, in the noise conditions 
SA— tSZ and SA+tSZ, the lowest string tensions where denois- 
ing is effective are greatly increased relative to the noise condition 
PA— IN. For SA— tSZ, the eye visibility threshold is slightly be- 
low the best experimental bound, while for SA+tSZ it is slightly 
above. These results are absolutely in the line of those obtained in 



the noise condition PA+IN, as the secondary anisotropics represent 
even stronger higher frequency noise. 



5.4 Comparison to PSM detectability threshold 

As our WDBD algorithm uses the PSM for preliminary localiza- 
tion of the string tension, it is a natural question to ask whether the 
overall denoising performance at low string tensions is limited by 
this preliminary PSM localization. We address this by defining and 
studying the detectability threshold for the PSM, which provides a 
measure of the minimum string tension where the PSM alone pro- 
vides robust detection of strings. 
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Figure 10. Posterior probability distributions (red solid curves) for the string tension as computed from a simulated signal observed at 1 arcminute resolution. 
The left and middle panels relate to the noise conditions SA— tSZ for string tensions p = 4.0 x 10"'^ and p = 1.0 x 10"'^ respectively. The right panel relates 
to the noise conditions SA+tSZ for a string tension p = 2.5xl0~^. The black dashed vertical lines represent the exact values of the string tension relative to 
the test string simulation. 





Figure 11. Signal-to-noise ratio (left panels) in decibels (dB), correlation coefficient (middle panels), and kurtosis (right panels) of the magnitude of gradient 
as functions of the string tension in logarithmic scaling in the range logio P ^ ["lO^ —05] and at 1 arcminute resolution. The top panels relate to the noise 
conditions SA— tSZ, while the bottom panels relate to the noise conditions SA+tSZ. The black dashed curves represent values before denoising, while the 
red solid curves and green dot-dashed curves represent values after WDBD and Wiener filtering respectively. The vertical lines on the curves represent the 
variability at one standard deviation of the estimated statistic across the 100 test simulations considered (these lines are not visible by eye where smaller than 
the width of the curves). The blue dotted vertical lines represent the eye visibility thresholds p = 1.0 x 10"'^ for the top panels and p = 2.5 x 10"'^ for the 
bottom panels. The blue dotted horizontal lines identify either the limit of zero signal-to-ratio, zero correlation coefficient, or the kurtosis of the magnitude of 
gradient of a pure string signal: ^cl^^l ~ 32. 



We firstly describe how the PSM detection threshold is com- 
puted, based on a hypothesis test for string detection. We may de- 
fine an estimation p of the string tension from the observed signal 
/ as the expectation value of the posterior probability distribution 
(I2TI) computed on the basis of the PSM: 

p = s[p(p|/)]. (28) 

For any possible string tension, the probability distribution function 
for ^may consequently be obtained from simulations. 

We identify the critical value po such that, for null string ten- 
sion, one has p{p ^ po) = a, for some suitable positive value a 



much smaller than unity. The test for the hypothesis of null string 
tension is then defined as follows. For estimated values p po, 
the hypothesis of null string tension may be rejected with a signif- 
icance level a. On the contrary for estimated values p < po, the 
hypothesis of null string tension may not be rejected. 

We define the detectability threshold p"^ such that, for a string 
tension p'^, one has p(p ^ po) — 1 — /3, for some other suitable 
positive value (3 much smaller than unity. Consequently, for string 
tensions larger than p* , the probability of rejecting a null string ten- 
sion on the basis of the hypothesis test defined is larger than 1 — (3. 
The value p^ is the smallest string tension that can be discrimi- 
nated from the hypothesis of null string tension for given values of 



© 2009 RAS, MNRAS 000,[TJ{T6l 



String signal denoising 15 



Noise condition 


PSM detectability 


Eye visibility 


PA-IN 


2.2x10-10 


6.3x10- 


■10 


PA+IN 


1.7x10-8 


2.5x10" 


-8 


SA-tSZ 


6.1x10-8 


1.0x10" 


-7 


SA+tSZ 


1.9x10-'^ 


2.5x10" 


-7 



Table 2. PSM detectability and eye visibility thresholds on the string ten- 
sion determined on the basis of the PSM for each of the noise conditions 
considered. All values are given with two significant figures. 

a and (3. It may thus be understood as a detectability threshold de- 
termined on the basis of the PSM. As our overall denoising method 
is using the PSM as a preliminary estimation of the string tension, 
identifies an effective lower bound on the string tension range 
where denoising could reasonably be expected to be effective. 

The PSM detectability thresholds in the various noise condi- 
tions considered are reported in Table |2] for a ^ [3 ^ 0.01. In 
all cases except SA+tSZ the PSM detectability thresholds are be- 
low the best experimental bound, while for SA+tSZ the PSM de- 
tectability threshold is around the best experimental bound. 

Secondly, we compare these thresholds with the eye visibility 
thresholds, which as noted previously also indicate the lower limit 
of string tensions where our quantitative measures show effective 
performance for the WDBD algorithm. For each of the noise condi- 
tions with significant high frequency content, i.e. PA+IN, SA— tSZ 
and SA+tSZ, the PSM detectability threshold p* is slightly below 
the eye visibility threshold. For the PA— IN case, this difference 
is larger, the PSM detectability threshold being about one third of 
the eye visibility threshold. This indicates that the PSM is able to 
more effectively exploit the high spatial frequency ranges where 
the string signal dominates the primary anisotropics. 

The discrepancy between these two thresholds shows that the 
detection problem alone can be solved with the PSM at slightly 
lower string tensions than the more difficult denoising problem. In- 
deed for values of the string tension between the two thresholds, 
denoising does not produce visible strings even though the PSM 
posterior probability distributions for p are distinctly peaked away 
from zero. It is one thing to estimate a single global parameter such 
as the string tension on the basis of a PSM, but quite another to 
explicitly reconstruct the string network itself. 

5.5 Algorithm robustness 

We comment here on the robustness of the WDBD algorithm rela- 
tive to both additional noise from foreground point sources and to 
the possible improvements in the definition of the denoising proce- 
dure itself. 

We have explicitly disregarded the problem of foreground 
emissions such as radio and infrared point sources. The discrim- 
ination of point sources from string loops imprinted in the CMB 
may appear to be a difficult task. However, the dipolar structure 
of the string loops represents an essential difference with point 
sources (Fraisse et al. 2008). In that context, the odd symmetry 
of the wavelets used in the WDBD algorithm (see Figure [3]) to 
match the string imprints is adequate both for long strings and 
string loops, and might help to discriminate between string loops 
and point sources. The algorithm was shown to be effective at de- 
tecting string loops, even at low tensions where long strings are 
not reconstructed anymore. However spurious point sources were 
also reconstructed at very small string tensions in the noise condi- 



tions PA— IN, in the absence of foreground point sources. A thor- 
ough analysis should be conducted in order to assess the real ro- 
bustness of the algorithm to discriminate between string loops and 
point sources, and to discuss necessary enhancements. 

Our approach explicitly assumes the statistical independence 
of coefficients of the wavelet decomposition, when conditioned 
on the string tension. However, significant correlations are present 
in the wavelet coefficients, and exploiting them should lead to 
improved denoising performance. Gaussian scale mixture (GSM) 
models may be considered which allow one to explicitly account 
for local correlations of the wavelet coefficients in the denoising 
process ([Andrews & Mallowsl [19741 : IPortilla et all 120031) . An en- 
hanced version of this model called the orientation- adapted Gaus- 
sian scale mixture (OAGSM) model relies on steerable wavelets 
in order to integ rate directional ity information in the local correla- 
tions ( Hammond & Simoncelli|[2008>) . A preliminary implementa- 
tion of the OAGSM model suggests that an enhancement relative to 
the WDBD algorithm may indeed be expected, albeit at significant 
computational cost. 

An improvement of the similarity of the shape of the filters 
to better match the typical string imprints may also be envisaged. 
Even though steerable wavelets can be very directional, their spatial 
support is not especially narrow. Filters with a more elongated sup- 
port such as curvelets ( Candes & Donoho 1999; Starck et al. 2 002h 
might be expected to provide better performance for the detection 
of long strings. Let us notice however that such filters would not be 
adequate anymore for string loops. Moreover a preliminary imple- 
mentation of this evolution provides no improvement relative to the 
WDBD algorithm for the detection of long strings. 

Finally, a discretization of the wavelet scales finer than the 
dyadic discretization used might provide an improved statistical 
model of the coefficients of the string signal at each spatial scale 
h. We did not consider this evolution here. 



6 CONCLUSION 

We have described a Bayesian framework for mapping the CMB 
signal induced by cosmic strings, based on a generalized Gaussian 
model capturing the sparse behaviour of the string signal in the 
steerable wavelet domain. This signal is buried in the standard pri- 
mary and secondary CMB anisotropics, which we model as Gaus- 
sian noise. For a fixed string tension we compute the Bayesian least 
squares estimator for each wavelet coefficient of the string signal. 
Our overall estimator is then formed as an average of these esti- 
mates for different string tensions, weighted by the posterior prob- 
ability of the string tension under a power spectral model. 

We have demonstrated the performance of our denoising al- 
gorithm through a series of numerical analyses at 1 arcminute res- 
olution consistent with upcoming experiments. The maps of the 
magnitude of the gradient of the denoised string signal produced 
by our algorithm were evaluated on the basis of three quantita- 
tive measures: the signal-to-noise ratio and correlation coefficient 
computed with respect to the known original string signal, and the 
kurtosis. In the idealized case of primary anisotropics without in- 
strumental noise, the strings can be identified for tensions down to 
p — 6.3xl0~^°, more than two orders of magnitude below the 
current experimental upper bound. With instrumental noise around 
1 /iK per pixel this lower bound is increased by more than one 
order of magnitude. The inclusion of secondary anisotropics fur- 
ther raises this bound to p — l.Ox 10~^ disregarding the thermal 
Sunyaev-Zel'dovich effect and to p = 2.5 x 10~^ including this 
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effect in the Rayleigh- Jeans limit. These values nonetheless remain 
slightly below or near the current experimental upper bound on the 
string tension, demonstrating that the proposed algorithm will be 
useful for analysis of upcoming high resolution data. 
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